LU factorization algorithm
The LU factorization algorithm was invented by the Polish mathematician Tadeusz Banachiewicz in 1938 (a professor of astronomy working at Jagiellonian University and AGH University in Krakow), which is described in [1].
Note that the Gaussian elimination algorithm generates an upper triangular matrix \( U \). It is also possible to generate a bottom triangular matrix \( L \) such that the original matrix of our problem \( A=LU \). What is the advantage of having such a decomposition? Suppose we want to solve two problems for two right-hand sides \( b_1, b_2 \). Let's also assume that we know the right side at the beginning \( b_1 \) and that the other right-hand side \( b_2 \) we will only know the other right-hand side b2 after solving the problem for the first right-hand side. Is there any way to use matrix factorization? \( A \) made for the previous right page \( b_1 \)? First we compute the LU factorization of the matrix \( A=LU \) then we have the equation \( LUx=b_1 \). We calculate first \( z \) by performing a forward substitution from equation \( Lz=b_1 \). Then we calculate x by solving the equation \( Ux=z \) by backward substitution.
Let us see how to modify the Gaussian elimination algorithm to generate both matrices \( L \) i \( U \). In order to obtain a matrix \( U \), we run the usual Gauss elimination algorithm (but without the right-hand side). We obtain Matrix \( L \) by remembering the numbers by which we multiply each line during elimination.
Initially
\( L= \begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ \end{bmatrix} \)
and
\( U= \begin{bmatrix} 0.05 & 0.10833 & 0.00833 & 0 & 0 \\ 0.10833 & 0.5 & 0.21666 & 0.00833 & 0 \\ 0.00833 & 0.21666 & 0.55 & 0.21666 & 0.00833 \\ 0 & 0.00833 & 0.21666 & 0.5 & 0.10833 \\ 0 & 0 & 0.00833 & 0.10833 & 0.05\\ \end{bmatrix} \)
During each multiplication, we update the matrix \( L \)
\( 1^{st} = 1^{st} \frac{1}{A_{1,1}}= 1^{st} \frac{1}{0.05}=1^{st}20 \)
\( U = \begin{bmatrix} {\color{blue}0.05 }*20 & 0.10833*20 & 0.00833*20 & 0*20 & 0*20 \\ 0.10833 & 0.5 & 0.21666 & 0.00833 & 0 \\ 0.00833 & 0.21666 & 0.55 & 0.21666 & 0.00833 \\ 0 & 0.00833 & 0.21666 & 0.5 & 0.10833 \\ 0 & 0 & 0.00833 & 0.10833 & 0.05 \\ \end{bmatrix} \)
In particular, we put the words from the diagonal of the original matrix on the diagonal.
\( L= \begin{bmatrix} {\color{blue}0.05} & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ {\color{blue}0.10833} & 0.5 & 0.21666 & 0.00833 & 0 \\ 0.00833 & 0.21666 & 0.55 & 0.21666 & 0.00833 \\ 0 & 0.00833 & 0.21666 & 0.5 & 0.10833 \\ 0 & 0 & 0.00833 & 0.10833 & 0.05\\ \end{bmatrix} \)
When subtracting lines, into the matrix \( U \) we insert the number used to scale the line \( 2^{nd} = 2^{nd} - 1^{st} *A_{2,1}= 2^{nd} - 1^{st} *{\color{blue}0.10833} \).
\( L= \begin{bmatrix} 0.05 & 0 & 0 & 0 & 0 \\ 0.10833 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ \end{bmatrix} \)
\( U \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0.10833-1*0.10833 & 0.5-2.16666*0.10833 & 0.21666-0.16666*0.10833 & 0.00833 & 0 \\ 0.00833 & 0.21666 & 0.55 & 0.21666 & 0.00833 \\ 0 & 0.00833 & 0.21666 & 0.5 & 0.10833 \\ 0 & 0 & 0.00833 & 0.10833 & 0.05\\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 0.26528 & 0.19860 & 0.00833 & 0 \\ {\color{blue}0.00833} & 0.21666 & 0.55 & 0.21666 & 0.00833 \\ 0 & 0.00833 & 0.21666 & 0.5 & 0.10833 \\ 0 & 0 & 0.00833 & 0.10833 & 0.05 \\ \end{bmatrix} \)
\( 3^{rd} = 3^{rd} - 1^{st} A_{3,1}= 3^{rd} - 1^{st} {\color{blue}0.00833} \)
\( L= \begin{bmatrix} 0.05 & 0 & 0 & 0 & 0 \\ 0.10833 & 1 & 0 & 0 & 0 \\ {\color{blue}0.00833} & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 0.26528 & 0.19860 & 0.00833 & 0 \\ 0.00833-1*0.00833 & 0.21666-2.16666*0.00833 & 0.55-0.16666*0.00833 & 0.21666 & 0.00833 \\ 0 & 0.00833 & 0.21666 & 0.5 & 0.10833 \\ 0 & 0 & 0.00833 & 0.10833 & 0.05 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 0.26528 & 0.19860 & 0.00833 & 0 \\ & 0.19861 & 0.54861 & 0.21666 & 0.00833 \\ 0 & 0.00833 & 0.21666 & 0.5 & 0.10833\\ 0 & 0 & 0.00833 & 0.10833 & 0.05\\ \end{bmatrix} \)
\( 2^{nd} = 2^{nd} \frac{1}{A_{2,2}}=2^{nd}\frac{1}{0.26528}=2^{nd}3.76963 \)
\( L= \begin{bmatrix} 0.05 & 0 & 0 & 0 & 0 \\ 0.10833 & {\color{blue}0.26528} & 0 & 0 & 0 \\ 0.00833 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & {\color{blue}0.26528}*3.76963 & 0.19860*3.76963 & 0.00833*3.76963 & 0 \\ 0 & 0.19861 & 0.54861 & 0.21666 & 0.00833 \\ 0 & 0.00833 & 0.21666 & 0.5 & 0.10833 \\ 0 & 0 & 0.00833 & 0.10833 & 0.05 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 1 & 0.74864 & 0.03140 & 0 \\ 0 & {\color{blue}0.19861} & 0.54861 & 0.21666 & 0.00833 \\ 0 & 0.00833 & 0.21666 & 0.5 & 0.10833\\ 0 & 0 & 0.00833 & 0.10833 & 0.05 \\ \end{bmatrix} \)
\( 3^{rd} = 3^{rd} - 2^{nd} A_{3,2}= 3^{rd} - 2^{st}{\color{blue}0.19861} \)
\( L= \begin{bmatrix} 0.05 & 0 & 0 & 0 & 0 \\ 0.10833 & 0.26528 & 0 & 0 & 0 \\ 0.00833 & {\color{blue}0.19861} & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 1 & 0.74864 & 0.03140 & 0 \\ 0 & 0.19861-1*0.19861 & 0.54861-0.74864*0.19861 & 0.21666-0.03140*0.19861 & 0.00833 \\ 0 & 0.00833 & 0.21666 & 0.5 & 0.10833 \\ 0 & 0 & 0.00833 & 0.10833 & 0.05 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 1 & 0.74864 & 0.03140 & 0 \\ 0 & 0 & 0.39992 & 0.21042 & 0.00833 \\ 0 & {\color{blue}0.00833} & 0.21666 & 0.5 & 0.10833 \\ 0 & 0 & 0.00833 & 0.10833 & 0.05 \\ \end{bmatrix} \)
\( 4^{th}=4^{th}-2^{nd}A_{4,2}=4^{th}-2^{nd}{\color{blue}0.00833} \)
\( L= \begin{bmatrix} 0.05 & 0 & 0 & 0 & 0 \\ 0.10833 & 0.26528 & 0 & 0 & 0 \\ 0.00833 & 0.19861 & 1 & 0 & 0 \\ 0 & {\color{blue}0.00833} & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 1 & 0.74864 & 0.03140 & 0 \\ 0 & 0 & 0.39992 & 0.21042 & 0.00833 \\ 0 & 0.00833-1*0.00833 & 0.21666-0.74864*0.00833 & 0.5-0.03140*0.00833 & 0.10833 \\ 0 & 0 & 0.00833 & 0.10833 & 0.05 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 1 & 0.74864 & 0.03140 & 0 \\ 0 & 0 & 0.39992 & 0.21042 & 0.00833 \\ 0 & 0 & 0.21042 & 0.499738 &0.10833 \\ 0 & 0 & 0.00833 & 0.10833 & 0.05 \\ \end{bmatrix} \)
\( 3^{rd}=3^{rd}\frac{1}{A_{3,3}}=3^{rd}\frac{1}{0.39992}=3^{rd}2.50050 \)
\( L=\begin{bmatrix} 0.05 & 0 & 0 & 0 & 0 \\ 0.10833 & 0.26528 & 0 & 0 & 0 \\ 0.00833 & 0.19861 & {\color{blue}0.39992} & 0 & 0 \\ 0 & 0.00833 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 1 & 0.74864 & 0.03140 & 0 \\ 0 & 0 & {\color{blue}0.39992}*2.50050 & 0.21042 & 0.00833*2.50050 \\ 0 & 0 & 0.21042 & 0.499738 & 0.10833 \\ 0 & 0 & 0.00833 & 0.10833 & 0.05 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 1 & 0.74864 & 0.03140 & 0 \\ 0 & 0 & 1 & 0.52615 & 0.02082 \\ 0 & 0 & {\color{blue}0.21042} & 0.499738 & 0.10833 \\ 0 & 0 & 0.00833 & 0.10833 & 0.05 \\ \end{bmatrix} \)
\( 4^{rd}=4^{rd}-3^{rd}A_{4,3}=4^{rd}-3^{rd}{\color{blue}0.21042} \)
\( L= \begin{bmatrix} 0.05 & 0 & 0 & 0 & 0 \\ 0.10833 & 0.26528 & 0 & 0 & 0 \\ 0.00833 & 0.19861 & 0.39992 & 0 & 0 \\ 0 & 0.00833 & {\color{blue}0.21042} & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 1 & 0.74864 & 0.03140 & 0 \\ 0 & 0 & 1 & 0.52615 & 0.02082 \\ 0 & 0 & 0.21042-1*0.21042 & 0.499738-0.52615*0.21042 & 0.10833-0.02082*0.21042 \\ 0 & 0 & 0.00833 & 0.10833 & 0.05 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 1 & 0.74864 & 0.03140 & 0 \\ 0 & 0 & 1 & 0.52615 & 0.02082 \\ 0 & 0 & 0 & 0.38902 & 0.10395 \\ 0 & 0 & {\color{blue}0.00833} & 0.10833 & 0.05 \\ \end{bmatrix} \)
\( 5^{th}=5^{th}-3^{rd}*A_{5,3}=5^{th}-3^{rd}{\color{blue}0.00833} \)
\( L= \begin{bmatrix} 0.05 & 0 & 0 & 0 & 0 \\ 0.10833 & 0.26528 & 0 & 0 & 0 \\ 0.00833 & 0.19861 & 0.39992 & 0 & 0 \\ 0 & 0.00833 & 0.21042 & 1 & 0 \\ 0 & 0 & {\color{blue}0.00833} & 0 & 1 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 1 & 0.74864 & 0.03140 & 0 \\ 0 & 0 & 1 & 0.52615 & 0.02082 \\ 0 & 0 & 0 & 0.38902 & 0.10395 \\ 0 & 0 & 0.00833-1*0.0.00833 & 0.10833-0.52615*0.00833 & 0.05-0.02082*0.00833 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 1 & 0.74864 & 0.03140 & 0 \\ 0 & 0 & 1 & 0.52615 & 0.02082 \\ 0 & 0 & 0 & 0.38902 & 0.10395 \\ 0 & 0 & 0 & 0.10394 & 0.04982 \\ \end{bmatrix} \)
\( 4^{th}=4^{th}\frac{1}{A_{4,4}}=4^{th}\frac{1}{0.38902}=4^{th}2.57056 \)
\( L= \begin{bmatrix} 0.05 & 0 & 0 & 0 & 0 \\ 0.10833 & 0.26528 & 0 & 0 & 0 \\ 0.00833 & 0.19861 & 0.39992 & 0 & 0 \\ 0 & 0.00833 & 0.21042 & {\color{blue}0.38902} & 0 \\ 0 & 0 & 0.00833 & 0 & 1 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 1 & 0.74864 & 0.03140 & 0 \\ 0 & 0 & 1 & 0.52615 & 0.02082 \\ 0 & 0 & 0 & {\color{blue}2.57056}*0.38902 & 2.57056*0.10395 \\ 0 & 0 & 0 & 0.10394 & 0,04982 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 1 & 0.74864 & 0.03140 & 0 \\ 0 & 0 & 1 & 0.52615 & 0.02082 \\ 0 & 0 & 0 & 1 & 0.26721 \\ 0 & 0 & 0 & {\color{blue}0.10394} & 0.04982 \\ \end{bmatrix} \)
\( 5^{th}=5^{th}-4^{th}*{\color{blue}0.10394} \)
\( L= \begin{bmatrix} 0.05 & 0 & 0 & 0 & 0 \\ 0.10833 & 0.26528 & 0 & 0 & 0 \\ 0.00833 & 0.19861 & 0.39992 & 0 & 0 \\ 0 & 0.00833 & 0.21042 & 0.38902 & 0 \\ 0 & 0 & 0.00833 & {\color{blue}0.10394} & 1 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 1 & 0.74864 & 0.03140 & 0 \\ 0 & 0 & 1 & 0.52615 & 0.02082 \\ 0 & 0 & 0 & 1 & 0.26721 \\ 0 & 0 & 0 & 0.010394-1*0.010394 & 0.04982-0.26721*0.010394 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 1 & 0.74864 & 0.03140 & 0 \\ 0 & 0 & 1 & 0.52615 & 0.02082 \\ 0 & 0 & 0 & 1 & 0.26721 \\ 0 & 0 & 0 & 0 & 0.04704 \\ \end{bmatrix} \)
\( 5^{th}=5^{th}-\frac{1}{A_{5,5}}=5^{th}\frac{1}{0,04704}=5^{th}*21.2585 \)
\( L= \begin{bmatrix} 0.05 & 0 & 0 & 0 & 0 \\ 0.10833 & 0.26528 & 0 & 0 & 0 \\ 0.00833 & 0.19861 & 0.39992 & 0 & 0 \\ 0 & 0.00833 & 0.21042 & 0.38902 & 0 \\ 0 & 0 & 0.00833 & 0.10394 & {\color{blue}0.047042} \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 1 & 0.74864 & 0.03140 & 0 \\ 0 & 0 & 1 & 0.52615 & 0.02082 \\ 0 & 0 & 0 & 1 & 0.26721 \\ 0 & 0 & 0 & 0\ & {\color{blue}0.047042}*21.2585 \\ \end{bmatrix} \)
As a result, we got the following matrices \( L \) i \( U \).
\( L * U = \begin{bmatrix} 0.05 & 0 & 0 & 0 & 0 \\ 0.10833 & 0.26528 & 0 & 0 & 0 \\ 0.00833 & 0.19861 & 0.39992 & 0 & 0 \\ 0 & 0.00833 & 0.21042 & 0.38902 & 0 \\ 0 & 0 & 0.00833 & 0.10394 & 0.047042 \\ \end{bmatrix} * \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 1 & 0.74864 & 0.03140 & 0 \\ 0 & 0 & 1 & 0.52615 & 0.02082 \\ 0 & 0 & 0 & 1 & 0.26721 \\ 0 & 0 & 0 & 0\ & 1 \\ \end{bmatrix} \)